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Via continuum mechanics [PRL 103,086401] with Random Phase Approximation (dRPA) screen- 
ing, we develop a numerically efficient general-geometry electronic exchange-correlation energy func- 
tional. It gives correct asymptotic power laws for dispersion interactions between insulators or met- 
als. As a numerical example we obtain the full binding energy curves e{D) for parallel metal slabs 
of small but finite thickness: at all separations D our e{D) agrees better with full dRPA correlation 
calculations than does the Local Density Approximation, while being much more efficient than full 
dRPA correlation. 
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An increasing body of workjll-Q has demonstrated 
that the correlation energy £'dRPA ^j^g direct Random- 
Phase Approximation (dRPA) is highly accurate for en- 
ergy differences in many and varied electronic systems, 
at least in cases where orbital self interaction is not an 
issue. dRPA binding properties for a wide variety of bulk 
materials [1] are typically more accurate than those from 
the local density approximation (LDA) , especially for dis- 
persion (van der Waals, vdW) bound systems !3]. For the 
vdW attractive potential, which is totally neglected in 
the LDA, the dRPA proves to be versatile, predicting 
unusual vdW coefficients @ and power lawsjg in agree- 
ment with quantum Monte Carlo results 0- 

^dRPA typically obtained in three steps: i) The 
bare response xo is obtained from occupied and un- 
occupied groundstate wavefunctions. This is typically 
the numerical bottleneck. Recent developments |8| at- 
tempt to bypass unoccupied states but can encounter 
problems for metallic systems, ii) The interacting re- 
sponse is calculated through the dRPA as xa(w) — 
Xo{(jj) + Xxo(u!)vx\{uj) where v is the Coulomb potential 
\r — r'\~^. iii) Finally the correlation energy is calculated 
via integration on the imaginary frequency axis through 
the Adiabatic Connection and Fluctuation Dissipation 
Theorem (ACFD) approach 
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where ^w) ~ — ■0^/^xo(<^)w^/^ is an Hermitian 
operator [sj. 

Other efficient van der Waals (vdW) functional [13, El 
give good results for many systems. However they rep- 
resent E^'^^ in an additive two-point approximation 
that is either obtained semi-empirically[l3j or derived [llj| 
by solving the dynamical screening problem ([T|) per- 
turbatively. As a result, these functionals miss non- 
pairwise-additive vdW energy contributions that can 
be substantial in highly polarizable, highly anisotropic 
systems jH, 0], including low-dimensional metals. Very 
large, anisotropic molecules and metallic and graphitic 
surface physics (e.g. binding of graphite on metal sur- 



faces) are two classes of systems where standard methods 
are inaccurate [l2| and dRPA is intractable. 

Here we solve equation ^ accurately thus avoiding the 
pairwise additive approximation, but we use the contin- 
uum mechanics of Tokatly, Tao, Gao and Vignale l^, lil 
to approximate xo in a numerically efficient manner. 
Their linearized continuum mechanics (CM) scheme 14 1 
uses the continuum fluid displacement u, which is related 
to the density perturbation by[lj, llSl j 



(2) 



For a small change to the Kohn-Sham (KS) potential 
V^{r,t) CM theory approximates u through the follow- 
ing hydro dynamic- like equation (from equations 3, 4 and 
14-16 of [3) 

~<^°„M^(r,t) + F^(r,t) 
n (r) 

where n°(r), = -71° {r)[d^,uV^^ {r)] and F°{r,t) de- 
pend on groundstate properties of the system. 

The force is defined in equation 14 of [l4|. Careful 



manipulation of equation 14 allows us to write it as = 
—Kf^^Ui^{r,t). Here K is a tensor, Hermitian [Kf^^ ~ 
KI ) operator defined by 



(4) 

(5) 
(6) 



It involves the electron density and groundstate ki- 
netic stress tensor = '^Y.i fAc>ti.'4^i{r)]*[duip.,{r)] ~ 
[d^„n (r)\ l^gj .^jjgj-g ^jjg gyjj^ is over occupied orbitals. 



In the absence of an external potential, equation ^ 
has time-periodic eigen-solutions defined by the hydro- 
dynamic eigen-equation 



(7) 



where N labels the sorted eigen- modes, UN{r) is related 
to an eigen-function of xo, > ^n-i is related to the 
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KS excitation energies (exactly in one-electron systems) 
and / drn°{r)u*^{r) ■ UM{r) = 6nm- 

By definition the tensor polarizability Xopi/(r, r'; w) is 
the time-periodic response of the /i cartesian component 
of the polarization —n'^{r)upf{r) to an external electric 
field in the v direction, while xo is the change in density 
n^(r) in response to a small change in the KS potential of 
form V^{r; r') = S{r—r'). They can be obtained through 
equations (H])-© and expansion in the eigen-solutions of 
([7]) provides the convenient forms (where in practise N 
is summed over the lowest A^Eig eigen-pairs) 



Xof^uir, r', ia) = ^ FNiicr)p%^{r)pNu(r') (8) 

N 

Xo{r,r',ta) = - ^ Fjv(*a)d^(r)djv(r') (9) 

N 

where F]\f{ia) = (17^ + cr^)^^, Pjv = n"{r)u]y{r) and 
djv(r) = [-V • {n°ir)uN{r)}] = -[V • pjvl- 

For efficient evaluation of the correlation energy, the 
following important relationships are derived which al- 
low us to evaluate Ec using integrals over one space 
variable only, reducing calculation time and storage 
requirements [1 7f . From ^ the projection of A [see ([T|)] 
in reciprocal space can be written in the separable form 
^J2NPN{icr)wl!{q)'WNiq') where 



WN 



(q) = - iq\Jv{q) ■ j dre"''''n°(r)Mjv(T") 



(10) 



or wn = y/v{q)dN{q) (here v{q) = Anq^^). Setting 
Wnai — J j0^wjM{q)w\j{q) allows us to define an 
A^Eig X -^Eig matrix ]B(i(T) with elements 
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with limjvEi,^ooTr^[G(B(iCT))] = Trg[G{A{ta))]'T^ for 
any analytic function G. 

Finally, defining the eigen- values oiM{i(T) to be (3i^{ia) 
we reduce the correlation energy ([T|) to the form 
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X ^Eo°g[i+/^''(*^)]-/3«(*^)}- (12) 



In practice we seem only to need a small number iVEig of 
eigen-solutions to converge correlation energies to a suf- 



ficiently small error (cx l/n% . ) within CM theory [18| 



This agrees with other observations (e.g. ref. |8|) that cal- 
culating Ec through a diagonalisation of xqv requires few 
eigenvalues for convergence. 

The most trying calculation in this functional method 
is evaluation of equation (O, as K is a spatially- 
dependent, differential operator. To overcome this prob- 
lem we use an auxiliary basis set B = {(/>j(T")}, of size 
A^Bas, which need not be mutually orthogonal but must 
be complete in the limit A^Bas — ^ oo. Choice of this basis 



is the only part of his scheme that differs for different 
geometries or systems: for example, plane waves for pe- 
riodic systems, gaussians for atoms and molecules. With 
a given basis set we expand our CM eigen-function ([7|) as 
UNfj.i'T') = '^Tv^t'^^i (*") which we substitute into equa- 
tion ([7]). This provides a set of SA^Bas x SA'^Bas coupled 
equations 

-n%N°ka%^ - {<f^V + K^k^..} (13) 

while Nji^{a-'j^^a'lf^) — 6nm sets the orthogonality. 

The non-operator terms in these equations 
are = J drn"{r)<j)*ir)(j3kir) and $0^^^ ^ 

- / dr[n°{r)dfj,„V^^{r)](p*{r)(l)k{r). Separating the 
final term into Kjk^^i, = / dr(p*{r)K^i,{r)(pk{r) = 



(T) 



— ^K^/l!,,, and using integration by parts gives 
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-HT;,[V0*]-[V0fc]}dr (14) 
: J n"[a,V0*] • [d^Wcl>k]dr (15) 



where all terms are functions of r and all derivatives can, 
ideally, be performed analytically on the basis functions. 

Surprisingly for a hydrodynamic-style approach, CM 
theory gives the exact bare responses Xq, xo to irrota- 
tional fields of one- and two(ti)-electron systems around 



their groundstate[13l |18|. This means that our corre- 
lation scheme will give the same results as dRPA for 
the asymptotic vdW interaction between two hydrogen 
or two helium atoms. 

To explore this further we follow [3| in expanding both 
the KS and CM response to 0(a;~^) leading to the fol- 
lowing identities [ist 



i=E^ 



jaN, 



(16) 



Here hjaN = '^^^^ where is a KS eigen- 

energy difference of an occupied orbital \j) and un- 
occupied orbital \a) with occupations fj and /a, and 
KjaN = / (j | 1^) ■ ujsidr is a mode-overlap matrix el- 
ement of the current (obtained via the current opera- 
tor J). Thus D,N > El — Eh for isolated systems and 
flNq > minfc(eLfc+q — en/e) for periodic systems where L 
labels the lowest unoccupied orbital or band and H labels 
the highest occupied. 

One implication of this is that a Kohn-Sham insulator 
will remain an insulator under CM, in the sense of finite 
responses dH)-® as ct ^ 0. Thus[lil CM theory obeys 
the well-known vdW laws for insulators with (e.g.) a 
— C4Z3^'* asymptotic binding for two thin layers. This 
is a very strong feature of the CM theory, not shared 
by common approximated ACFD theories [III [l9| where 
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explicit cutoffs have to be imposed in the tails in order 
to suppress metallic-like response. 

In the opposite limit of a homogeneous electron gas 
(HEG), CM is analytically soluble, agrees with the true 
Xo for q <Si kp, uj 3> vpQ^ and in particular has a "metal- 
lic" infinite polarizability, Xq — > oo as g and cj — > 0. 
Electron-gas-like (metallic) systems nevertheless pose a 
difficult test for CM theory because the single-particle- 
like excitations occurring for u < vpq, (and thus not 
accurately desribed by CM) , can make significant contri- 
butions to the RPA correlations, mainly at short spatial 
range (large wavelength). 

This inaccuracy can be improved in metallic systems 
by employing range-separation (RS) such that the short- 
range physics is treated by a local scheme. This makes no 
contribution to vdW asymptotic physics. A well-studied 
RS scheme is described in It involves choosing a 

gns and splitting up the Coulomb potential, with a long- 
range component v^'^^^'^{r) — erf(gRsr)r~^ , equivalent to 
replacing dTO]) by w^'''-'(g) = W7v(9)e-9'/(s«^s). We la- 

IrCM 

bel the corresponding correlation energy i?c This 
has the additional benefit of accelerating convergence. 

For Xo to be reliably approximated by continuum me- 
chanics without a separate treatment of the low frequen- 
cies we must choose g^s to be substantially less than kp. 
Here we use g^s = 0.25rJ^ — O.lSkp where is a global 
measure of the inter-electron distance. For the jellium 
slab problems studied below we simply choose corre- 
sponding to the background charge density of each slab, 
though more general prescriptions exist. The remaining 
correlation must be included from local approximations 
so that 



CM 



where gj^'^''^''^^' is the correlation energy per electron of 
the HEG with a short-ranged interaction, taken from ^21]. 

Ideally we must also implement a range-separation for 
exchange, but this proves numerically difficult for the 
slab geometries we investigate. We instead use the ra- 
tio of the long-range exchange to total exchange of an 
HEG Ax « l.lgRs^s/ v^l + (l-l'ZRS^s)^ as a prefactor 
for the exact exchange (EXX) E^^^ = -i J drdr'\r - 
r'l^^l ^„ /„?/',* ('")V'n('"')P ^'^id make up the remainder 
with the LDA. Combining this with (fT7|) gives i?xc = 
A.£;Exx + (1 _ Ax)E^^^ + E^^'^-^' . 

As a numerical test of our proposed functional we 
choose the difficult case of two thin metal slabs described 
in Refs. [l^ and [2^ This system is defined by three 
parameters only: the width of the slabs s, the inner 
surface-surface distance D and the positive background 
charge electron density p — 3/(47rr^). The total number 
of electrons per unit area is Ng — 2sp ~ rt+(z)dz. 

We test our method on slab pairs with s = 3ao, = 
1.25ao and s = 5ao, Ts = 2.07ao which have been studied 
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FIG. 1. e(D) graph for = 1.25, s = 3. RPA data from [23 
Inset data shows the vdW dominated region. 



in Ref. [231 and Refs. [19| andl22| respectively. Especially in 
the first case the LDA and dRPA give significantly differ- 
ent energy curves. We consider the cleavage energy per 
electron £(£>) = e{D)~e{oo) = [Eo{D) ~ Eo{oo)]/Ns as a 
function of D. Slabs with rs < 4 have a defined binding 
length Dq where the force is zero. Thus a binding energy 
£(, = [^(ZJo)! and an elastic modulus C'zz — dnDf^iDo) can 
also be defined. 

In Figure[T]we plot e{D) versus D for = 1.25, s — 3. 
Our method matches the RPA closely for this system. 
Binding properties for both studied systems are tabu- 
lated in Table U and show that the = 2.07, s = 5 
system is less well-predicted but still much better than 
the LDA. If instead we set gRs — oo the results be- 
come much worse for both cases. For widely separated 
slabs {D :s> s) the CM theory correctly and analytically 
describes coupled two-dimensional plasmons and hence 
correctly predicts the known asymptotic dRPA formf^ 
e{D > s) w -0.012562v^(£'+.s)"'^/^ With s = 12.8ao 
and Vs = 2ao . . . 6ao we calculate C5 /2 numerically within 
8% of the theory. By contrast most other efficient vdW 
functionals would predict an incorrect power law expo- 
nent in this limit with e{D) « — C4Z?~*. 

In our CM calculations we use auxiliary basis func- 
tions (/)fcq||(r) = 6fc(z)e~"^ii'''«fT8^ at iV,,, qy points. All 
calculations are quite efficient with the slowest step be- 

). Convergence 



ing evaluation of Wnm at 0{Nq^^ ^q^^b. 



is reached with iVsa 



42, N„ 



55, N„ = 250 and 





CM LDA dRPA 


CM LDA dRPA 




rs = 1.25, s = 3 


Ts = 2.07, s = 5 


Do 

Czz 


3.33 3.38 3.32t 
0.74 0.53 0.79t 
0.51 0.45 0.55t 


1.57 1.56 1.62±0.1§ 
1.78 1.72 1.85±0.1§ 
1.31 1.38 1.32±0.1§ 



TABLE I. Groundstate properties of two slab systems under 
different approximations. Energies are in mHa/e~ and dis- 
tance are in Bohr radii. J from Ref. [23 . § is guessed from 
Refs. [T9I and taking into account estimated error bars. 
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N^ig < 60. Our dRPA calculation takes approximately 
eight times longer than the groundstate LDA calculation. 
Test runs of full RPA calculations for these systems took 
hours, compared to minutes for our functional, consistent 
with A^qii — 55. We also note that a 10% variation in q^g 
made only a 1% change to e^. 

While results for our test systems are not perfect, they 
show closer agreement with the dRPA than the LDA both 
in the binding region and for larger D, with a marked im- 
provement in speed over full dRPA. The vdW dispersive 
physics is treated accurately and shows excellent agree- 
ment with the dRPA in contrast to other methods. The 
current prescription has a wide scope for refinement both 
empirically through adjustment of gRs and the exchange 
functional and by introducing better physics, most ob- 
viously through improved (semi-local) treatment of low- 
frequency behaviour which will reduce dependence on the 
range separation. 

Furthermore preliminary tests suggest that metals are 
a worst-case for CM theory - i.e. that range separation 
will be much less needed for bound and insulating sys- 
tems. 

In summary, we have derived and developed an ef- 
ficient general-geometry functional with correct long- 
range correlation energy. Its ability to predict correctly 
the vdW physics of metallic and insulating systems is a 
distinct advantage over other efficient vdW functionals. 
It is currently being implemented for periodic systems, 
which should enable meaningful energy calculations for 
(e.g.) vdW bonded nanosystems such as metallic nan- 
otube arrays or graphene on metals. These systems 
require non-pair-additive high-level computations (e.g. 
dRPA) with a large unit cell, which is beyond present 
computational power. 

The authors would like to thank I. Tokatly, J. Jung, 
A. Savin, J. Angyan, , and G. Vignale for fruitful discus- 
sions. 
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EQUIVALENCE OF TRACES 

The calculation of Ec via W relies on an equivalence of 
certain traces. Here we demonstrate that 

TVq[log(l + A)-A] =TrAr[log(I + B) - B] (1) 

where A = -vkji^f^yk^ Bnm = \J Fn {io)FM Wnm 
and Wnm = / ■0^WN{q)wlf{q). 

Let us project A onto Fourier space giving 

Aiq,q')^Y.^¥UMwl,iq)wMiq') (2) 

NM 

where F^va/ — FN{ia)SNM and we subsequently omit a 
for notational brevity. We can then evaluate 
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or more generally 



AM{q,q') ^ - Y,^*NiQ)wM{q')W' 
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(3) 
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when M{q,q') = Ewm ^"7v(g)wAf (gOMwM. 
Observation of A^ suggests that 

AP{q,q')^Y.'^*Niq)wM{q'WWr-^¥]NM. 

NM 

Premultiplying by A validates this assumption since 
AAP{q,q') =J2wNiQ)wMiq'WW{¥W)P-^¥ 
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~-AP+\q,q') 
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and thus ([5]) is the correct form by induction on p. 
Finally we can take traces (and use permutations under 



\NN 



a trace) to show 

Tr,[i]-E / j^^'M^NiqWh 
= YWnnW]nn 

N 

=TrAr[FW] = TrAr[B] (7) 
Tr,[AP]=J2 I ^w*Aq>M{qWWY-'¥UM 
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^Ti-n[{¥W)p] = TrATpf] 



(8) 



where 



VFWVF and B 



NM 



y/ Fn {ia)FM (icr) Wn m ■ Finally it follows from Taylor 
expansion of an analytic function G{A) — gpA^ that 



Tr,[G(i)]=^gpTr,[i'']=Tr 



TV 



(9) 



Setting G{x) — log(l + x) — x gives the equivalence. 



EIGENVALUE CONVERGENCE 

In the main work we mention the rapid convergence 
with respect to the CM eigen-pairs. While this is some- 
what system specific some general considerations make it 
likely that fewer eigen-pairs than Kohn-Sham (KS) tran- 
sition modes will be required to obtain the same conver- 
gence in most systems. 

For any given KS problem we can expand solutions in 
a finite basis set (eg. a real space grid, planewaves, Gaus- 
sians etc.) of size iVeas which can be made as large as we 
like. The KS equations thus have Nbus solutions of which 
A^Occ niay be considered occupied and Nbus — Nqcc are 
unoccupied. Since involves a sum over occupied and 
unoccupied states, calculations involving x have a lead- 
ing 0(iVocc(^Bas - Nocc)) ~ 0(7Vocc^Bas) if wc assumc 
A^Occ < Nb^. 

By contrast the tensor CM equations have N-^ig = 
SA^Bas solutions. Projection of xo involves a sum over 
all solutions and takes 0(3iVBas)- Combining this with 
the sum rules [see (16) of main text] means that the CM 
eigcn-frequcncies will be distributed more sparsely than 
the KS-transitions (except in one and two electron 
systems). 
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In the specific case of a periodic system this has a very 
notable effect on integration over the Brillioun zone. Here 
a natural basis set for both the LDA and CM is e*('J+<^) '' 
where q lies within the Brillioun zone and G is a recipro- 
cal lattice vector. It is thus sufhcient to project xo onto 
q, G and G' . Following the standard notation we use k 
rather than q for the groundstate sampling. 

In a full KS-dRPA calculation, projection of xo is 
0{Nocc bandA^Band-^^fc^c) whcrc A^Qcc band IS the num- 
ber of occupied bands and A^Band is the total number of 
unoccupied bands. is the number of points sampled 
in the Brillioun zone and must be counted twice: once for 
each projection onto q and once to cover all transitions 
from fc to fc -I- q (see e.g. section III of ref. [H for a more 
comprehensive discussion). By contrast the same projec- 
tion in CM is 0{NEigNkN^) as all A; to A; -I- q transitions 
are treated collectively in a single uj^g eigenmode. 

This relationship can also be observed if we solve for x 
by direct perturbation Q instead of the diagonalization 
method used here. In a full KS calculation we need to 
perturb each of the Aqcc orbitals to obtain x- In the CM 
we need only to solve for the three-dimensional quantity 
u. In the periodic bulk case A'occ = N^Nocc band- 

With a similar energy cutoff for both the KS and CM 
solutions so that A^Eig ~ A^Band, a planewave-based eval- 
uation of the CM-dRPA will run at least 0{Nk) times 
quicker than a full dRPA calculation. In an efficient 
dRPA code in an insulating system, A'fc is typically found 
to be converged with between 200 and 2000 points^l] di- 
vided by the number of symmetries. 

Under the dRPA, the truncation error caused by using 
only a finite number of transitions, is 0(1/cj^^j^) where 
'^a^jm — ^a„^ ~ (where am is an occupied KS state 
and jm is unoccupied) is the largest transition frequency 
included in the sum. The equivalent error for the CM is 
0(1/J77v„) where i^N^ is the transition frequency of the 
highest included mode. 

We must note, however, that there are a variety of 
ways to calculate Ec- For example we can bypass the 
projection of XOi as we do in the method presented in 
the manuscript, by using W. For periodic bulks this is an 
0{NkN^^Nc) calculation under the CM with 0{N^^^) 
storage and 0{N^^^) diagonalisation. In the full KS 
it would be OiNlN^^^Ni^^Nc) with OiN^N^^^NiJ 
storage and 0(A^;^A^occ-^Bas) diagonalisation which is in- 
feasible for most bulk systems. 

As a general rule, projection of xo onto space will be 
more efficient than using W for a true KS bulk or peri- 
odic system. Calculation of W will be quicker for the CM 
in almost all geometries and will be quicker than evalu- 
ation of W for the full KS system in all but the smallest 
systems. 

In summary, for systems with more than a handful 
of occupied orbitals, the CM is all but guaranteed to 
converge faster than the full KS under the same method 
for calculating Ec. Furthermore it may offer the potential 



to use alternate methods that are faster still (such as 
using W for correlation energies) which may be infeasible 
or much slower in the full KS system. 

EXACTNESS IN ONE-ELECTRON SYSTEMS 

Equations (41), (45-46) of Ref. i provide a proof 
that the continuum mechanics approach is exact for 
one-electron (or two-electron with equal spin up and 
down) systems. This follows from the equivalence of 
CM to Madelung dynamics which are equivalent to the 
Schrodinger equation for one-orbital systems in irrota- 
tional fields jj]. We will provide a direct proof from the 
Schrodinger equation elsewhere. 

NON-CONTRIBUTING GAUGE MODES 

It is possible to find solutions mjv of the CM equations 
where V • hPum — 0- For example, in the 2DEG systems 
examined modes of form un — g_Le"^" ^ where q^ = 
q|l X 2: have this property. 

For such modes (labelled A^*) the weights of the sum 
rules [equation (16) of main text] become hjaN* = 0. 
While a full proof is beyond the scope of this sup- 
plement it follows from the exactness of the CM re- 
sponse to 0(lli~*) and the near-completeness over all oc- 
cupied/unoccupied pairs ja of {j\J\a) for vectors with 
non-zero gradients. Thus J2ja ^jaN' — and = 
J2ja^jaN'0J^a = 0- cousequcncc of this is that 

any CM eigen-mode of a KS-insulating system which has 
Hn' — must also have V • hF'um' = 0. 

Such modes will either be supressed by the boundary 
conditions of the problem, or not contribute to the corre- 
lation energy of the system and may thus be discarded. 
More precisely they do not contribute to the linear re- 
sponse Xo at all and contribute only a shift of gauge 
to the tensor response Xq. The former identity can be 
seen directly as only V • itPun appears in the expression 
[equation (9) of main text] for xo while the latter comes 
indirectly from integration by parts over the Coulomb 
tensor. 

INTERACTING INSULATOR 

We show in the main work that a KS insulating system 
(at a bare response level) will remain insulating under the 
CM. More formally this means that both the KS system 
and the CM approximation thereto, obey 

\Xog{G,G'■^a)\<^^^^^ (10) 

where Y and W are undetermined but non-zero con- 
stants. In the CM case this follows from Ppf{r) = 
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n°{r)uN{r) < oo, PAr(r oo) = and n^q > and 
the asymptote linio— j-oo |Xoq| = Cja^ . 

Using the and matrix representation of the 
RPA-interacting response Xa it is possible to show 



H NM 



where ^^f^^ is an eigenvalue of the matrix D^q = 
diag[17^q] + AWq. Provided IS)\q does not have any zero 
eigenvalues (this would be a strongly correlated metalli- 
sation and is extremely rare) it is obvious that 



\Xxq{G,G';ia)\ < 



for some finite Y' and W' . 



Y' 



(12) 



THIN INSULATORS HAVE D * POWER LAWS 

If we have two well-separated electronic systems at 
a distance D with no electronic overlap between them, 
we can write the dRPA dispersion energy in Lifshitz-like 
form: 

rvdW/ 



da 



T{a) 

T{a) =Tr log(l - XiaVabXibVba) 



(13) 
(14) 



where Xia/b is the interacting response of system A/B 
in isolation defined in the dRPA as 



XlA/B —XOA/B + XOA/B^XIA/B- 



(15) 



Vab — Vba is the Coulomb potential between electrons 
in different systems only and can be considered a function 
of rAjfB and D where r^/^ is a position in system A/B. 
Equations ([T3|)-([T4| follow from analysis of Feynmann 
ring diagrams (see [5| for a similar analysis) or from the 
ACFD where the intra-system Coulomb interaction v and 
Vab are switched on separately. 

Let us define a system composed of two periodic slabs 
localised in z (ie. where we can define a length s such that 
XiA/Bq{z, z') is negligible for \z\ > s/2 or \z'\ > s/2) and 
centered such that J zdznP^^^{z) — 0. In such a system 
it can be shown[6] that, in the limit _D 3> s and D — >■ oo, 
(fT4)) is equivalent to 



T{a) 
Here 



BZ 



dqii log(l-Xi^q||Xisq||g|e-29|i^) . (16) 



X 



/dz/d2:'(q|| - iz) ■ 
XiA/Sqii (0) -2, 0, z'; id) • (^11 + iz) where we project 
the tensor response XiA/Bq{G\\, z,G'^^, z';ia) in recipro- 
cal space in the xy-plane and real space in z. 



In an insulating system we have shown that |Xiq (icr)! < 



Y 



and thus |Xiq(z(T)| < 



Y' 



where Y and thus 



Y' is finite. Since W is finite there is a well-defined up- 
per bound for insulating systems and we need not worry 
about singularities in XiA/Bq we can set x = q\\D and 
make a series expansion in Z?~^ such that 



27ra;da;log[l - XiAo^iBoj^e ^] 



^XiAoX^Bo + OiD-"^). 



Inserting this into (jl6p we find 



da 



XiAoiia)XiBo{ia-) 



(17) 

(18) 
(19) 



Ergo two insulating slabs have a I?"'' power law in CM 
theory as in full dRPA. This conclusion applies only for 
insulators: for thin metals the a integration diverges. 



TWO-SLAB GEOMETRY 



Let us define our two-slab metal problem to have a 



background charge n+(z) ~ p[H{-, 



z-L 



ff(f - 



\z + L\)] where L = and H{x) = IVx > 0, other- 

wise. This defines two jellium slabs of width s, surface-to- 
surface distance D and backround charge per unit area 
p — 3/(47rr^). The total number of electrons per unit 
area is set to Ng — 2s p — n+(z). 

The partial isotropy means V^^{r) = V^^{z) and the 
KS wavefunctions take the form 



^nk, (r) =p„(z)e-''=ii-''ii 



(20) 



where / dzp* (z)pm(z) = (27r)^^(5,„„ . The KS energies 
are e„fc|| — eno + ^|fc||P with occupation /„ — 2max(ei? — 
e„o, 0). The density and kinetic pressure tensor are thus 



T{z) =t°^z)[x <g) X + y <g) ij] +t^'%z)z(giz (22) 



(21) 
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|p„(z)|2 and t°-(z) = 



where t°\\{z) = E„ 
E„/J9,p„(z)p-ia,,n0(z). 

For the present slab problem we choose auxiliary basis 
functions of the form (/)fcq|| (r) — 6fe(z)e~*''ii '''n where 6fc(z) 
is either or tanh(fctz)e*"~ where kt is a parameter 

chosen to optimise convergence and n is an integer. The 
restiction to integer n makes this basis set incomplete 
but inclusion of non-integer n does not alter results. 

We then set 



iNq. 



k 



0fcqiiM[aL(9||)^ + a^||(9||)g||] (23) 
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(the q±= q\\ x i: term does not contribute to the corre- 
lation energy). Thus the eigen-equations are 

^lin ) =Kjk\\ II (911 (9ii ) 

+ K,k\iAq\\VNAq\\) (24) 

+ K,kz\\iq\\)a%\iiq\\) (25) 
which must be solved for each q||. Normalisation gives 

Here N^k = J dzn°{z)b*{z)bk{z) and = 
/ dzn°(2;)[9^zy^^(2;)]6*(2;)6i;(z) are independent of 
The components of K((7|| ) take the form 

-Kjkzz =3-Fjfc[t°", 1, 1] + i-^jfe[n", 2, 2] 

+ g2 (^J-^.,[tO-,0,0] + ij-,fe[n", 1, 1]^ (26) 

+ gf (^3J-jfe[t°ll,0,0] + ij-jfe[n",l,l]) (27) 
=(-«g||) +t°ll,0,l] + ij-jfcK,l,2]^ 

(28) 



and 



+ -^-^jfeK,o,i]. 



where we use the shorthand J-jk[f,a,b] 

J fiz)[d,ab*{z)][dMz)]dz. 

Finally in this basis 



X ydze'«^^[g||Ujvq|||| +g2UArq||^], (29) 
WNM{q\,) = - I ir^w^ iqz)wMq„ (Qz)- (30) 



da f 27rg||dq|| 



27r J (27r)2 



Tr[L(B((j||,za))] (31) 



Bnm ^VfNi<y)fMia)WNMiq\\) (32) 



where L(x) — log(l + a;) — a;. 

In our calculations we use approximately 500-1000 reg- 
ularly distributed z points for quadrature (with the num- 
ber depending on system size). We also use approx- 
imately 500 Qz points on a Gauss-Hermite grid (due 
to the Range-Separation term e~'^ /(29rs) calculate 
WNMiq\\)- This is more than sufficient to represent the 
chosen basis functions in either space. 

To correctly integrate over frequency we require a grid 
that accurately deals with functions of form a/(fe^ + u^) 
where b ranges from very small to large. Choosing a 
regular grid for a <^ 1 and using a Clenshaw-Curtis grid 
for larger a seems to work well for these problems. 
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